#include <assert.h>
#include<mpi.h>
#include<stdlib.h>
#include <iostream>
using namespace std;

typedef unsigned long int INDEX;

int init() { return 0;}
int accum(int t, int d) { return t+d;}
int combine(int left, int right) { return left+right;}
int scan(int t, int not_used) { return t;}

int main(int argc, char* argv[]){
  int np;
  int rank;
  int ierr;
  MPI_Status status;
  INDEX INPUT_SIZE;
  INDEX size;

  int myTally;
  int lTally;
  int myTallyt;
  int pTally;
  int stride = 1;
  int* data;
  int* result;

  // MPI init
  MPI_Init(&argc, &argv);
  MPI_Comm_size(MPI_COMM_WORLD, &np);
  MPI_Comm_rank(MPI_COMM_WORLD, &rank);

  // Load file & init
  if(rank==0){
    assert(argc > 1);
    FILE* fin=fopen(argv[1],"r");
    fscanf(fin,"%ld", &INPUT_SIZE);

    data = (int*) new int [INPUT_SIZE];
    result = (int*) new int [INPUT_SIZE];

    for(INDEX i=0;i< INPUT_SIZE;i++){
      fscanf(fin,"%d", &data[i]);
    }
    fclose(fin);
  }

  // Broadcast input size
  if (MPI_Bcast(&INPUT_SIZE, 1, MPI_UNSIGNED_LONG, 0, MPI_COMM_WORLD) != MPI_SUCCESS) {
    cout<<"MPI_Bcast failed"<<endl;
    fflush(stdout);
    MPI_Abort(MPI_COMM_WORLD, -1);
  }

  // Allocate memory for local process
  size = INPUT_SIZE/np;
  int* ldata = (int*) new int [size];
  int* lresult = (int*) new int [size];

  // Distribute data for local process
  ierr = MPI_Scatter(data, size, MPI_INT,
      ldata, size, MPI_INT,
      0, MPI_COMM_WORLD);

  // Local process
  myTally = init();
  for(INDEX i=0; i<size; i++){
    myTally = accum(myTally, ldata[i]);
  }

  // Going up
  while(stride < np) {
    //lTally[rank+stride] = myTally[rank];
    //myTally[rank] = combine(myTally[rank], myTally[rank+stride]);
    MPI_Barrier(MPI_COMM_WORLD);
    if (rank%(2*stride) == 0) {
      MPI_Send(&myTally, 1, MPI_INT, rank+stride, 0, MPI_COMM_WORLD);
      MPI_Recv(&myTallyt, 1, MPI_INT, rank+stride, 0, MPI_COMM_WORLD, &status);
      myTally = combine(myTally, myTallyt);
    }
    if ((rank-stride)%(2*stride) == 0) {
      MPI_Recv(&lTally, 1, MPI_INT, rank-stride, 0, MPI_COMM_WORLD, &status);
      MPI_Send(&myTally, 1, MPI_INT, rank-stride, 0, MPI_COMM_WORLD);
    }
    stride = 2*stride;
  }

  if (rank == 0) {
    myTally = init();
  }

  // Propogate down
  stride = np/2;
  while (stride >= 1) {
    MPI_Barrier(MPI_COMM_WORLD);
    //myTally[rank+stride] = combine(myTally[threadid], lTally[rank+stride]);
    if (rank%(2*stride) == 0) {
      MPI_Send(&myTally, 1, MPI_INT, rank+stride, 0, MPI_COMM_WORLD);
    }
    if ((rank-stride)%(2*stride) == 0) {
      MPI_Recv(&myTallyt, 1, MPI_INT, rank-stride, 0, MPI_COMM_WORLD, &status);
      myTally = combine(myTallyt, lTally);
    }
    stride = stride/2;
  }
  MPI_Barrier(MPI_COMM_WORLD);
  pTally = myTally;
  for(INDEX i=0; i<size; i++) {
    pTally=accum(pTally, ldata[i]);
    lresult[i] = scan(pTally, ldata[i]);
  }
  if(rank==0){
    for(INDEX i=0; i<size; i++) {
      result[i] = lresult[i];
    }
    for(int i=1; i<np;i++){
      MPI_Recv(result+i*size, size, MPI_INT, i, 0, MPI_COMM_WORLD, &status);
    }
  }else{
    MPI_Send(lresult, size, MPI_INT, 0, 0, MPI_COMM_WORLD);
  }

  // Clean up
  delete [] ldata;
  delete [] lresult;
  if(rank==0){
    for(INDEX i=0;i< INPUT_SIZE;i++){
      printf("%d ",result[i]);
    }
    printf("\n");
  }
  MPI_Finalize();
  return 0;
}
